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ABSTRACT 


We review the difficulties of the classical fission and fragmentation hypotheses 
for the formation of binary and multiple stars. A crucial missing ingredient in 
previous theoretical studies is the inclusion of dynamically important levels of 
magnetic fields. As a minimal model for a candidate presursor to the formation 
of binary and multiple stars, we therefore formulate and solve the problem of 
the equilibria of isopedically magnetized, singular isothermal disks, without the 
assumption of axial symmetry. Considerable analytical progress can be made 
if we restrict our attention to models that are scale-free, i.e., that have surface 
densities that vary inversely with distance zu from the rotation axis of the system. 
In agreement with earlier analysis by Syer and Tremaine, we find that lopsided 
(A/ = 1) configurations exist at any dimensionless rotation rate, including zero. 
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Multiple-lobed (M = 2, 3, 4, ...) configurations bifurcate from an underlying 
axisymmetric sequence at progressively higher dimensionless rates of rotation, 
but such nonaxisymmetric sequences always terminate in shockwaves before they 
have a chance to fission into M = 2, 3, 4, ... separate bodies. On the basis of our 
experience in this paper, and the preceding Paper I, we advance the hypothesis 
that binary and multiple star-formation from smooth (i.e., not highly turbulent) 
starting states that are supercritical but in unstable mechanical balance requires 
the rapid (i.e., dynamical) loss of magnetic flux at some stage of the ensuing 
gravitational collapse. 


Subject headings: Hydrodynamics, Magnetohydrodynamics. Molecular Clouds, 
Stars: Binaries, Stars: Formation 


1. Introduction: Figures of Equilibrium and Binary Star Formation 

1.1. The Fission Hypothesis 

The fission hypothesis for binary star formation evolved from Newton’s calculation in 
the seventeenth century for the shape of a rotating Earth. Newton imagined an ingenious 
experiment boring holes to the center of our planet and filling them with water to show 
that the Earth is flatter at the poles than at the equator. This conclusion embroiled him in 
controversy with Cassini, who claimed on the basis of astronomical measurements that the 
Earth is prolate rather than oblate. (See Todhunter 1873 for a more detailed description, 
in particular, for an account of Maupertuis’s expedition to Lapland that settled the debate 
empirically in favor of Newton.) 

Newton's analysis assumed that the gravitational field of a homogeneous spherical Earth 
is undistorted by its slow rotation, with the centrifugal effects taken into account only in the 
fluid equilibrium . The general analytic expression describing the self-consistent eccentricity 
e = \J\ — £3/ £[ of an equilibrium spheroid of uniform density p with principal axes £3 < £2 = 
that rotates with constant angular velocity Q. was given by Maclaurin in 1742: 

02 2C1— p 2 ) 1/2 fi 

0 = ^ g~ p = — jr ~' ( 3 ~ 2e2 ) sin "‘ e - - e2 ) • (!) 

In the following year, Simpson (more widely known in connection with his “rule”) noticed 
that the Maclaurin spheroids can exist only if the rotational parameter /? < 0.449331. For 
P less than this critical value, two solutions exist, one more flattened than the other. At 


3 = 0, these two solutions correspond to a sphere (e = 0, most easily imagined in the limit 
^ 0 with p finite) and a razor thin disk (e = 1, most easily imagined in the limit p — > oo 

with surface density E = f pdz and Q finite). 

Ninety one years later, Jacobi (1834) became intrigued by the existence of two entirely 
separate equilibria at low 0. He was particularly impressed by the fact that the less-obvious 
disk-like solution cannot be accessed from the spheroidal solution by means of a linear 
perturbation analysis. The presence of two unrelated solutions suggested to him that others 
may also exist. Jacobi relaxed the requirement of axisymmetry and showed that uniformly 
rotating, self-gravitating, liquid, masses can also assume triaxial equilibrium figures in which 
the principal axes £ y , and £ 3 have unequal values. 

Meyer (in 1842) discovered that the Jacobian sequence of triaxial ellipsoids branch- 
es from the Maclaurin spheroids when the latter’s eccentricity reaches e = 0.81267 (0 = 
0.37423). At that point, the figure axes i\ and i 2 of the Jacobian ellipsoids become equal, 
and Jacobian sequence merges into the Maclaurin sequence. If a Maclaurin spheroid is 
allowed to dissipate energy and contract homologously to higher density while conserving 
angular momentum, it will become triaxial before e can exceed 0.81267. In other words, 
the Maclaurin spheroids are secularly unstable with respect to viscous forces and bifurcation 
into Jacobian ellipsoids 1 . 

In 1885, Poincare found that the Jacobian sequence bifurcates into further classes of 
equilibrium that have lop-sided shapes. The first bifurcation sequence corresponds to a series 
of egg-shaped figures that become pear-shaped, and occurs when 0 = 0.28403. Poincare 
envisioned the slow evolution of a contracting spheroid in which the contraction time scale is 
much longer than the internal viscous timescale so that uniform rotation can be maintained. 
Such an object was imagined to progress along the Maclaurin sequence as it spins up. Upon 
reaching 3 = 0.3r423, it would lose its axial symmetry and become a Jacobian ellipsoid. 
Poincare then conjectured that further secular evolution to 3 = 0.28403 and beyond would 
lead to bifurcation into the pear-shaped sequence of figures, which, in the face of additional 
increases in the density and rotation rate, would eventually fission into a parent body and 
a satellite, such as the Earth and its Moon. The same sequence of events was invoked by 
G.H. Darwin (1906), the son of the naturalist, to account for the origin of binary stars (see 
also Darwin’s 1909 review). 


Consult Chandrasekhar (1969) for an account of the dynamical instability of Maclaurin spheroids a- 
gainst transformation into Riemann ellipsoids that contain internal circulation. He also analyzed the secular 
instability of rotating ellipsoids against transformation by gravitational radiation into Dedekind ellipsoids 
whose figure axes remain fixed in space. 



Liapounoff (1905), Jeans (1916), and Cartan (1928), however, discovered that the Ja- 
cobi sequence of ellipsoids becomes dynamically unstable at exactly the point (/? = 0.28403) 
where Poincare s pear-shaped figures first appear. The inevitable appearance of dynamical 
instabilities renders the fission hypothesis problematical, in part because of the mathematical 
difficulties associated with describing three-dimensional nonlinear hydrodynamical evolution. 
A more fundamental difficulty arises from uniformly rotating gaseous equilibrium configura- 
tions with realistic degrees of central condensation (for example, gaseous polytropes) reaching 
equatorial breakup prior to bifurcation into triaxial configurations (James 1964). Further- 
more, if, as likely, internal viscous timescales exceed the contraction timescale, a polytropic 
configuration will develop differential rotation. As clarified by Ostriker & Mark (1968), and 
Ostriker & Bodenheimer (1973), contracting differentially-rotating polytropes become bar- 
unstable before reaching equatorial breakup. Therefore, a realistic modern descendant of the 
fission hypothesis would amount to the conjecture that an unstable barred figure fragments 
into two or more pieces. This hypothesis foundered when definitive numerical simulations by 
Durisen et al. (1986) demonstrated that the emergent bar drives spiral waves that transport 
angular momentum outward and mass inward, in the process stabilizing the configuration 
against fission. Astronomically, this result is consistent with the observation that bars in 
flattened galaxies drive outer spiral structures, and do not spin off additional galaxies. 


1.2. The Fragmentation Hypothesis 

An alternative theory for the formation of binary stars can be traced back to Jeans 
(1902), who specified the minimum mass, Mj oc G _3 / 2 a 3 p -1 / 2 for an object of isothermal 
sound speed a and mean density p, to collapse under its self-gravity in the presence of op- 
posing gradients of gas pressure (see also Ebert 1955, and Bonnor 1955) . Hoyle (1953) 
considered a large cloud with mass M ~ Mj initially. As it collapses, with a held constant 
(because radiative losses under optically thin conditions tend to keep cosmic gases isother- 
mal) but p increasing, the cloud progressively contains additional Jeans-mass subunits, which 
might collapse individually onto their own centers of attraction. Adjacent collapsing subfrag- 
ments could then conceivably wind up as binary stars. A stability analysis by Hunter (1962) 
of homogeneously collapsing, pressure-free spheres seemed to support the Hoyle conjecture. 
However, Lavzer (1964) argued that because the overall collapse and the growth of perturba- 
tions proceed with the same powers of time, individual subunits may have insufficient time 
to condense into independent entities before the entire cloud disappeared into the singularity 
of Hunter s background state (the analog of the big crunch in a closed-universe cosmology). 

A further difficulty with the fragmentation hypothesis arises because self-gravitating 


systems that are initially close to hydrostatic equilibrium (or have only one Jeans mass) are 
necessarily centrally condensed. Numerical calculations by Larson (1969) indicated that such 
centrally condensed masses would collapse highly non-homologously. In the case of a singular 
isothermal sphere - which has a density distribution p = a 2 /2nGr 2 and which contains one 
Jeans mass at each radius r - Shu (1977) showed that collapse proceeds in a self-similar 
manner, from inside-out . Past the moment t — 0 when collapse is initiated, a rarefaction 
wave moves outward at the speed of sound a into the hydrostatic envelope of the cloud. At 
an y given time t > 0, roughly half of the disturbed material is infalling, and half has been 
incorporated into a tiny hydrostatic central protostar approximated as a mass point. At no 
time in the process does any subvolume excluding the center contain more than one Jeans 
mass. Shu (1977) conjectured that such solutions are unlikely to fragment, a conclusion 
verified by Tohline (1982) to apply more generally to a wide variety of centrally-condensed 
collapses. 

If such a collapsing cloud is imbued with angular momentum, a structure containing 
a star/disk/infalling-envelope naturally develops (Terebey, Shu k Cassen 1984). Numerical 
work by Boss (1993) removing the assumption of axial symmetry indicates that rotating 
collapse flows with radial density profiles as centrally concentrated as p a r' 2 also avoid 
fragmentation on the way down. The fragmentation hypothesis is therefore restricted either 
to cases of the collapse of less centrally condensed clouds (e.g. Burkert, Bate k Bodenheimer 

1997). or else to cases of breakup into multiple gravitating bodies after a disk has already 
formed. 

Although the issues of gravitational instabilities and fragmentation within disks are 
still active areas of investigation, calculations by Laughlin k Bodenheimer (1994), which 
specifically followed the nonaxisymmetric evolution of disks arising from the collapse of 
rotating r 2 clouds, did not find disk fragmentation (see also Tomley et al. 1994; Pickett, 
et al. 1998). Rather, as the disks arising from the collapse flow become gravitationally 
unstable, they develop spiral structures which elicit an inward flux of mass and an outward 
flux of angular momentum that proves sufficiently efficient as to stabilize the disk against 
fragmentation (see also Laughlin, Korchagin k Adams 1998). 

Boss (1993) has conjectured that isolated molecular cloud cores with density laws as 
steep as p a r~ 2 will inevitably lead to the formation of single stars accompanied by planets 
rather than binary stars. Since most stars in the Galaxy are members of multiple systems, he 
concludes that collapsing cloud cores must generally arise from configurations less steep than 
p <x r This point of view is supported by Ward-Thompson et al. (1994), who claim that 
observed prestellar molecular cloud cores always have substantial central portions that are 
flat, p ^ const, rather than continue along the power law, p oc r“ 2 , that characterizes their 



outer regions. It should be noted, however, that such configurations are in fact consistent 
with the predictions of theoretical calculations of molecular-cloud core-evolution by ambipo- 
lar diffusion (Nakano 1979, Lizano & Shu 1989, Basu & Mouschovias 1994), which show that 
nearly pure power-laws, p oc r -2 , arise only for a single instant in time, the pivotal state (Li k 
Shu 1996), just before the onset of protostar formation by dynamical infall. Moreover, more 
recent analyses of the millimeter- and submillimeter-wave dust-emission profiles by Evans et 
al. (2000) and Zucconi et al. (2000) that take into account the drop in dust temperature 
(but perhaps, not the gas temperature) in the central regions of externally irradiated dark 
clouds show that the portion of the density profile of prestellar cloud cores that is flat (p « 
const), if present at all, is considerably smaller than originally estimated by Ward-Thompson 
et al. (1994). 

One can also note that while the Taurus molecular-cloud region represents the classic 
case of isolated star formation (Myers k Benson 1983), it contains, if anything, more than 
its cosmic share of binaries (Ghez, Neugebauer k Matthews 1993; Leinert et al. 1993; 
Mathieu 1994; Simon et al. 1995; Brandner et al. 1996). Moreover, when observed by 
radio-interferometric techniques, Taurus contains many cloud cores that are well fit by p oc 
r -2 envelopes, yet each star-forming core typically contains multiple young stellar objects 
(Looney. Mundv k Welch 1997). 

Recent high-resolution simulations of the fragmentation problem carried out with a- 
daptive-mesh techniques (Truelove et al. 1998) indicate that many of the previous hydro- 
dvnamical simulations claiming successful fragmentation with density laws less steep than 
p oc r~ 2 contained serious errors. Indeed, as long as the starting conditions are smooth and 
close to being in mechanical equilibrium (i.e. , start with only one Jeans mass), gravitational 
collapses seem vn general not to produce fragmentation. The emphasis on the sole fault lying 
with the law p oc r~ 2 is therefore misplaced. Something else is needed. Klein et al. (2000) 
identify the missing ingredient as cloud turbulence; our opinion is that magnetic fields may 
be equally or even more important. 


1.3. The Effect of Magnetic Fields 

It is a proposition universally acknowledged that on scales larger than small dense cores, 
magnetic fields are more important than thermal pressure (but perhaps not turbulence) in the 
support of molecular clouds against their self-gravitation (see the review of Shu, Adams, & 
Lizano 1987). Mestel has long emphasized that the presence of dynamically significant levels 
of magnetic fields changes the fragmentation problem completely (Mestel k Spitzer 1956; 
Mestel 1965a, b; Mestel 1985). Associated with the flux <I> frozen into a cloud (or any piece 


of a cloud) is a magnetic critical mass: 


McrW = 


2wG 1 


( 2 ) 


Subcritical clouds with masses M less than A/ cr have magnetic (tension) forces that are 
generally larger than and in opposition to self-graviration (e.g., Shu k Li 1997) and cannot 
be induced to collapse by any increase of the external pressure. Supercritical clouds with 
M > M cr do have the analog of the Jeans mass - or. more properly, the Bonnor-Ebert mass 
- definable for them, but unless they are highly supercritical, M » M cr , they do not easily 
fragment upon gravitational contraction. The reason is that if M ~ M cr for the cloud as a 
whole, then any piece of it is likely to be subcritical since the attached mass of the piece scales 
as its volume, whereas the attached flux scales as its cross-sectional area. Indeed, the piece 
remains subcritical for any amount of contraction oi the system, as long as the assumption 
of field freezing applies. An exception holds if the ' loud is highly flattened, in which case 
the enclosed mass and enclosed flux of smaller pieces both scale as the cross-sectional area. 
This observation led Mestel (1965, 1985) to specula in that isothermal supercritical clouds, 
upon contraction into highly flattened objects, could and would gravitationally fragment. 
The present paper casts doubt on this speculation (a , when the original cloud begins from a 
state of mechanical equilibrium, and (b) when magnetic; flux is conserved by the contracting 
cloud (see also Shu Li 1997). 


Zeeman observations of numerous regions (see the summary by Crutcher 1999) indicate 
that molecular clouds are, at best, only marginally •>U|)ercritical. The result may be easily 
justified after the fact as a selection bias (Shu et al. 1999). Highly supercritical clouds have 
evidently long ago collapsed into stars; they are not found in the Galaxy today. Highly 
subcritical clouds are not self-gravitating regions; they must be held in bv external pressure 
(or by converging fluid motions); thus, they do not constitute the star-forming molecular- 
clouds that are candidates for the Zeeman measurements summarized by Crutcher (1999). 
The clouds (and cloud cores) of interest for star formation today are, by this line of reasoning, 
marginally supercritical almost by default. 

The above comments motivate our interest in re-examining the entire question of binary- 
star formation by the fission and fragmentation mechanisms, but including the all-important 
dynamical effects of magnetic fields and the empirically well-founded assumption that pre- 
collapse cloud cores have radial density profiles that, in first approximation, can be taken 
as p ex r . Li & Shu (1996: see also Baureis, Ebert k Schmitz 1989) have shown that 
the general, axisymmetric, magnetized equilibria representing such pivotal states assume the 
form of singular isothermal toroids (SITs): p(r,9) oc r *R(9) in spherical polar coordinates 
(r,9,(p), where R(9) = 0 for 9 = 0 and n (i.e., the density vanishes along the magnetic 
poles). YYe regard these equilibria as the isothermal (rather than incompressible) analogs 
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of Maclaurin spheroids, but with the flattening produced by magnetic fields rather than by 
rotation. In the limit of vanishing magnetic support, SITs become singular isothermal spheres 
(SISs). In the limit where magnetic support is infinitely more important than isothermal 
gas pressure, SITs become singular isothermal disks (SIDs), with p(w,z) = E (vj)5(z) in 
cylindrical coordinates (zu ,<p,z), where S(z) is the Dirac delta function, and the surface 
density E(co) oc £7~ l . 

In a fashion analogous to the SIS (Shu 1977), the gravitational collapses of SITs have 
elegant self-similar properties (Allen k Shu 2000). But it should be clear that the formation 
of binary and multiple stars could never result from any calculation that imposes a priori an 
assumption of axial symmetry. In this regard, we would do well to remember the warning of 
Jacobi in 1834: 

"One would make a grave mistake if one supposed that the axisymmetric spheroids of revo- 
lution are the only admissible figures of equilibrium.” 

Motivated by the insights of those who have preceded us, we therefore start the cam- 
paign to understand binary and multiple star-formation by considering in this paper the 
nonaxisymmetric equilibria of self-gravitating, magnetized, differentially-rotating, complete- 
ly flattened SIDs, with critical or supercritical ratios of mass-to-flux in units of (27 r(7 1/ ' 2 ) -1 . 

A = 2t (3) 

with A > 1 (see Li k Shu 1996, Shu k Li 1997). Keeping A fixed, i.e., under the assumption 
of field freezing, we shall find that such sequences of non-axisymmetric SIDs bifurcate from 
their axisymmetric counterparts at the analog of the dimensionless squared rotation rate /? 
(which we denote in our problem as D 2 ) given by the linearized stability analysis of Paper I 
(Shu et al. 2000; see also Sver k Tremaine 1996). Although some of these (Dedekind-like) 
sequences produce buds that look as if they might separate into two or more bodies, we 
find that, before the separation can be completed (by secular evolution?), the sequences 
terminate in shockwaves that transport angular momentum outward and mass inward in 
such a fashion as to prevent fission. 

In a future study, we shall follow the gravitational collapse of some of these non- 
axisymmetric pivotal SIDs. The linearized stability analysis and nonlinear simulations of 
Paper I suggests that the collapse of gravitationally unstable axisymmetric SIDs lead to 
configurations that are stable to further collapse but dynamically unstable to an infinity of 
nonaxisymmetric spiral modes that again transport angular momentum outward and mass 
inward in such a fashion as to prevent disk fragmentation. We suspect the same fate awaits 
the collapse of pivotal SIDs that are non-axisymmetric to begin with, as long as we continue 
with the assumption of field freezing. Thus, we shall speculate that rapid (i.e., dynamical 


rather than quasi-static) flux loss during some stage of the star formation process is an es- 
sential ingredient to the process of gravitational fragmentation to form binary and multiple 
stars from present-day molecular clouds. 

The rest of this paper is organized as follows. In §2 we derive the general equations 
governing the equilibrium of magnetized, scale-free, non-axisymmetric. self-gravitating SIDs 
with uniform velocity fields. In §3 we show that for SIDs with no internal motions the equa- 
tions of the problem can be solved analytically. For the more general case, in §4 we present 
an analytical treatment of the slightly nonlinear regime, when deviations from axisymmetry 
are small, valid for arbitrary values of the internal velocity field. In §5 we describe a numer- 
ical scheme to compute non-axisymmetric SIDs for arbitrary values of the parameters of the 
problem. Finally, in §6 we summarize the implications of our findings for a viable theory of 
binary and multiple star-formation from the gravitational collapse of supercritical molecular 
cloud cores that start out in a pivotal state of unstable mechanical equilibrium. 


2. Magnetized Singular Isothermal Disks 


The governing equations of our problem are given in Paper I (see also Shu & Li 1997). 
Thej are the usual gas dynamical equations for a completely flattened disk confined to the 
plane 2 — 0 except for two modifications introduced by the presence of magnetic fields that 
thread vertically through the disk, and that fan out above and below it without returning 
back to the disk. 


First, magnetic tension reduces the effective gravitational constant by a multiplicative 
factor e < 1, where 


with the dimensionless mass-to-flux ratio A > 1 taken to be a constant both spatially (the 
isopedic assumption) and temporally (the field-freezing assumption). Second, the gas pres- 
sure is augmented by the presence of magnetic pressure; this increases the square of the 
effective sound speed by a multiplicative factor 0 > 1, where we follow Paper I in adopting 


A 2 + 

A 2 + 1 ’ 


( 5 ) 



2.1. Equations for Steady Flow 


Consider the time-independent equation of continuity in 2D: 

V • (Eu) = 0. 

This equation can be trivially satisified by adopting a streamfunction defined by 

Eu = V x (tfe,), 

which written in cylindrical polar coordinates reads 

_ 1 9 $ 1 3 $ 

otE dip U<f E dzj 

Notice that u • V'E = 0, so curves of constant # describe streamlines. 


( 6 ) 

( 7 ) 

( 8 ) 


The momentum equation along streamlines can be replaced by Bernoulli’s theorem: 


^|u| 2 -I- 0"H(E) + eV = 5(\f), 


( 9 ) 


where B('i) is the Bernoulli function and ' H(E ) is the specific enthalpy associated with a 
barotropic equation of state (EOS) for the gas alone: 

-s dU. dE 




(10) 


In equation (10) the vertically intgrated pressure n is assumed to be a function of surface 
density E alone. For an isothermal EOS. we have n = a 2 E with a 2 = const, so that 
= a 2 In E plus an arbitrary additive constant that we are free to specify for calculational 
convenience. 


In terms of the variables introduced above, the vector momentum equation can now be 
written 

(Vxulxu + B'fljV^fl. (11) 


Expressed in component form, this equation gives the additional independent relation for 
momentum balance across streamlines: 


i r 5 




dzo 


y E dw ) + dip \ E dp> J 


= E &(V). 


(12 


Notice that the LHS is the z-component of — V x u: thus, E B' is the local vorticity contained 
in the flow (proportional to Oort’s B constant). The above set of equations is closed by the 
addition of Poisson’s equation: 


V(t*7, ip) = -G 



E(r, 0) rdr 

r 2 + zo' 2 — 2rrz;cos(0 — v 3 )] 1 ^ 2 


( 13 ) 


2.2. Scale-Free Isothermal Solutions 


For aligned SIDs, we look for solutions of the form, 


W(E) = a 2 lim In ( ) , 

R-+oo \ K / 

(14) 

E(zz,tp) = —S(<p), 

(15) 


where the constant K with dimension of g cm -1 and the dimensionless function S(<p) are to 
be determined. In equation (14) and in everything that follows, the limit operation R — » oo 
is to be taken after differentiation of variables like H and V in the equations of motion have 
occurred. We have taken advantage of the fact that additive constants in variables like R, 
V, and ^ do not enter the physical equations of motion to introduce a temporary artificial 
radial scale R so that we need not take logarithms of dimensional quantities. Putting the 
freedom to scale E entirely into K , we are free to normalize the function S(ip) such that 

^ f S(ip)dtp= 1. (16) 

Substitution of equation (15) into equation (13) yields 

V(«7, ip) = -GK Jim l S (iP)diP ( R/ \ a --- dX — , (17) 

R-*ooJ Jo [1 — X 2 — 2xCOs(^r — 7/.’)] 1 / 2 ' 7 

where x = r/zz. The inner integral can be evaluated by elementary techniques and gives 

ln{(R/zu) — cos (<p — xj;) + [1 + (R/zz) 2 — 2 {R/zj) cos(<^ - ^)j 1/2 } — ln[l — cos(<^> - r/>)]. (18) 

The argument of the first logarithm equals 

{R/zz) jl - ( zz/R ) cos (<p - ij') + [l - 2 (ct/ R) cos(^ - ip) + ( zzjR ) 2 J 1/2 j , (19) 

which can be expanded for large R as 

{2R/zj){\ - (zz/R) cos (cp - ip) + . . .]. (20) 

Thus, the inner integral in equation (17) equals 

ln(2 R/zz) + ln[l - (zz/R) cos (if - ip) + . . .} - ln[l - cos(</? - ip)]. (21) 

For large R , the middle logarithm goes to zero, and the substitution of the above result then 
yields 

V(zz,ip) = 27 tGI\ lim \\n(zz/'2R) + V(ip )\ , 


( 22 ) 



vv 


here 


V(<p) = ^ j> S(ip) In [ 1 - cos(t p - <p)] dip. 


We further look for solutions of the form, 


'&{vo,ip) = Q x/2 aK lim [-D \n(w/2R) + W{ip)), 

B(V) = - 0 1/2 a 

K 


(23) 


(24) 


(25) 


where D and B are dimensionless constants whose values are yet to be specified. In what 
follows, it is convenient to define the dimensionless radial mass flux as 


U(tp) = w'(<p). 


(26) 


which we will regard as an ODE for the angular part of the streamfunction XV (<p) if we know 
U(<p). An integration of equation (26) over a complete cycle shows that the mass flow across 
a full circle must vanish, 

jlJ(<p)d<p = 0, (27) 

since W(<p) is a periodic function of c p. In order for equation (27) to hold nontriviallv, U (ip) 
must possess both positive and negative values; thus, it must pass through zero at least once 
in the range (— tt, +tt). We define our angular coordinate so that U(<p) is zero at (p = 0: 

U( 0) = 0. (28) 

This convention results in U(<p) being an odd function of ip. 

Substitution of the expression for E and $ into equation (8) now yields the identifica- 


tions: 


U T . 


= e x/2 a 


u = e 1/2 q - 

s(p) ' * S(ip) 


(29) 


In other words, apart from the compression and decompression factor S(<p) as fluid elements 
flow in azimuth in a nonaxisymmetric disk, the dimensionless function U(tp) is the generator 
for radial motions and the dimensionless constant D is the generator for angular motions. 

The substitution of equations (10), (15), (23), (24), (25), and (29) into equation (9) now 
yields I< from the the radial part of the equality, 

0a 2 


K 

2TT6G 

whereas the angular part of the equality gives 

1 


(1 + DB), 


(30) 


2 S' 2 


(l J 2 + D 2 ) + (1 + DB)V + In 5 = -BW. 


(31) 


Similarly, equation (12) leads to the requirement, 

D d (U\ 

S + dp (5) 5 ' 

Since the combination U/S must be a periodic function of y ?, we may integrate equation (32) 
over a complete cycle and obtain the further constraint: 


27T J S(tp) 


(33) 


Finally, differentiating equation (31) with respect to tp and using equation (32) we obtain 


(S 2 - D 2 )S' + DUS + (1 + DB)S 3 V' = 0. (34) 

Equation (34) possesses critical points at 5(y?) = D, where u^, becomes equal to the magne- 
tosonic speed (see eq. 29). 


Equations (23), (32), (33) and (34) are the fundamental set of integro-differential equa- 
tions governing the problem. They have to be solved in the interval tp = [0, 2tt] for the three 
unknown functions S(p), V(p), U(p) and the unknown constant B. The constant D itself 
is freely specifiable. Notice that the arbitrarily introduced radial scale R enters nowhere in 
the final equations. 


Notice also that equation (32) implies that radial motions arise only in response to a 
local imbalance of forces — gravitational, pressure, and inertial — across streamlines, even 
though equation (33) requires such forces to be balanced on average over a circle. Moreover, 
the governing equations (23), (32), and (34) require S(y?) and V{p) to be symmetric with 
respect to p = ir when U{tp) is chosen to be antisymmetric. In other words, £(y?) and V{p) 
are cosine series in tp when 0(p) is developed as a sine series. Consequently, the choice of 
the zero of the angular coordinate is not unique: for a configuration with a basic .\/-fold 
symmetry, where M is a positive integer, the condition U{p) = 0 is satisfied in the interval 
[0, 2 tt] at p = kir/M with k = 0, 1, 2, ... , 2 M, and different choices of the x-axis correspond 
to rotations of the equilibrium configuration by multiples of 7 r /M. 


2.3. Fourier Decomposition for Poisson Integral 

\\ hen we have departures from axial symmetry, the difference form of the kernel and the 
periodic nature of the wanted solutions makes equation (23) suitable for solution by Fourier 
series. Since S(y?) and l (p) are periodic functions of p, they are expandable as the series: 

OO 

S(p) = S 0 + S m cos(mp). 

m— l 


( 35 ) 



( 36 ) 


u 


oo 

V(if) = Vo + 51 KnCOs(mc^), 

m— 1 

where we have isolated the axisymmetric terms So and V 0 . The coefficients S m and V m are 
real for all m > 1. In writing S(<p) and V((p) as pure cosine series, we have made use of our 
freedom to orient one of the principal figure axes of aligned equilibria along the z-axis. 

Since equation (23) gives V(tp) as a convolution of S(<p) and ln(l - cos<p), substitution 
of equations (35) and (36) into equation (23) and application of the convolution theorem for 
Fourier cosine transforms result in the identification, 

v;„ = L m S m , (37) 

where 

If. . . . , f —In 2 if m = 0 , ooX 

L m = -fHl-co S <p)cos(m V ,)d V =l ( _ l/lm \ ifW£1 (38) 

as shown in the Appendix. Therefore, the normalization condition (16) implies 

So = 1, VJ, = — In 2, (39) 

and equations (35) and (36) can be written as 

OO 

S(<^) = 1 — 5Z m ^m cos(nup). (40) 

m= 1 

oo 

V(ip) = -ln2+ 51 V m cos(rrup). (41) 

m=l 

If we are given {V' m }“ =1 , then we know S(y?) and V(<p). Unfortunately, local knowledge 
of either S or V at <p does not determine the value of the other at the same p. The relationship 

is local in Fourier space, so only global knowledge of S(<p) gives global knowledge of V'(v>); 

i.e., S(</?) is a functional, and not a function, of V(ip). 

Notice that in general, if a set V m of Fourier coefficients corresponds to a solution, then 
the set ( — 1 ) m V m corresponds to the same configuration rotated bv an angle n, as discussed 
at the end of §2.2. 


3. Static Equilibria 

For static equilibria, U = D = 0 and equation (34) reduces to 


5' + SV' = 0, 


(42) 


which has the barometric solution: S(p) = Ae~ v( ‘ p \ where A is a constant that can be 
adjusted to satisfy the normalization condition (16). Substitution of S(p) = .4e _l ^ into 
Poisson’s integral (eq. 23), or alternatively its Fourier decomposition (eq. 40 and 41), then 
constrains the solution for V '(</?). Remarkably, this system of nonlinear functional relations 
has an analytical solution, where iso-surface-density contours are ellipses of eccentricity e, 


%>) 


s/T^ 


(43) 


lie cos ip ’ 

with 0 < e < 1, the ± sign representing our freedom to rotate the equilibrium configuration 
by an angle 7 r, as anticipated at the end of §(2.2). Let us consider the case with the minus 
sign (the proof for the case with the plus sign is completely analogous). 


The set of Fourier coefficients corresponding to equation (43) is 

2\A 


V m = — <fi S(p) cos (nup) dp> = - - — - — — f 

7r m J 7r m Jo 1 — 


cos (rrup) 2 

dp = 

e cos p m 


T - y/l = 


(44) 

where we have used formula (3.613.1) of Gradshteyn & Ryzhik (1965) to evaluate the integral. 
^,From equation (41) we obtain the potential 


V(p) = -ln2-2 £ 

m=l 


1 - VI- e 2 


cos(m<p) 


in 


= In 


1 — e cos ^ 

1 + v 7 ! - e 2 ' 


(45) 


where we have evaluated the sum of the series by using formula (1.448.2) of Gradshteyn 
& Ryzhik (1965). By direct substitution, therefore, we see that S and V are related by 
the barometric relationship implied bv equation (42): S = Ae~ v , where .4 = \/l — e 2 /(l + 
VT^). (QED) 


The static axisymmetric solution (a magnetized but nonrotating disk with surface den- 
sity E = K/rv) is trivially recovered setting e = 0; Li k Shu (1997) give the time-dependent 
self-similar gravitational collapse of this special case. In the other extreme, for e = 1 the 
potential becomes 

V(p) = ln(l - cos<p), (46) 

and the corresponding set of Fourier coefficients, V m = —2/m, substituted into equation 
(40), gives the familiar Fourier expansion of the Dirac 5-function, 

OG 

S(p) = 1 + 2 cos(m<p) = 2 nS(p). (47) 

m= 1 

Thus, the equilibrium configuration degenerates into a semi-infinite filament with uniform 
linear mass density p = 2i tK . For values of e between these two extremes, both iso-surface- 
density contours and equipotentials are confocal ellipses of eccentricity e. Figure 1 shows 
some examples of static SIDs for different values of the eccentricity e. 



3.1. A Specific Example: the Molecular Cloud Core L1544 


As an amusing sideshow, Figure 2 shows an overlay of one of our eccentrically displaced 
static models projected onto a map of thermal dust emission at 1.3 mm obtained by Ward- 
Thompson, Motte, & Andre (1999) for the prestellar molecular cloud core L1544. Apart 
from relatively minor fluctuations due to the cloud turbulence, the solid curves depicting the 
iso-surface-density contours of the theoretical model match well both the observed shapes 
and grey-scale of the dust isophotes. 


Zeeman measurements of the magnetic-field component parallel to our line of sight 
toward LI 544 have been made by Crutcher & Troland (2000), who obtain B\\ = 11 ± 2 fiG. 
For a highly flattened disk, which is reflection symmetric about the plane 2 = 0, integration 
along the line of sight yields cancelling contributions of B~ and B v to By. The 2 -component 
of the magnetic field of our model core is given by 


We may now calculate the average value of E within a radius R as 

0a 2 A 2 (A 2 -f- 3) a 2 


(E) = <f dip [ Etz? dw = 

7 T tl J JO 


7 xeGR (A 4 — 1) -nGR' 


(48) 


(49) 


where we have made use of equations (4), (5), (15), (16) and (30). Therefore, the average 
value of B z within a radius R is 


( B : ) 


2t rG 1 / 2 
A 


<E> 


A(A 2 + 3) 2a 2 
(A 4 - 1) G X ! 2 R' 


(50) 


Notice the pleasant result that the above formulae do not involve e. 


Since we model LI 544 as a thin disk with elliptical iso-surface-density contours, its 
orientation in space is defined by three angles, two specifying the orientation of the disk 
plane, the third giving the position of the elliptical contours in this plane. We fix the first 
angle by assuming for simplicity that the major axis of the elliptical contours lies in the 
plane of the sky. The second angle i is the inclination of the minor axis with respect to 
the plane of the sky ( i = 0 for a face-on disk) and can be adjusted to fit the observations. 
The third angle, specifying the ellipse’s orientation in the disk plane, is given as 38° north 
through east by Ward-Thompson et al. (2000). 

We choose the eccentricity e and inclination i by the following procedure. From Figure 
2, we can estimate that a typical dust contour has a ratio of distances closest and farthest 
from the core center given in a model of nested confocal ellipses by 

~ 0.30 => 


1 + e 


e ~ 0.54. 


(51) 



Similarly, we may estimate that these ellipses have an apparent minor-to-major axis-ratio of 

(1 - e 2 )^ 2 cos i « 0.54 => cos i m 0.64. (52) 

The resulting ellipses for three iso-surface-density contours, spaced in a geometric progression 
1:2:4, are shown as solid curves in Figure 2. 

Determination of cosz allows us to compute an expected (Z?||) = {B z ) cost. Similarly, 
we obtain the expected hydrogen column density by multiplying (£) by (cosz) -1 for a slant 
path through an inclined sheet and by 0.7 for the mass fraction of H nuclei of mass m^: 
Nh — 0.7(E)/(m// cos z). 

The sound speed for the 10 K gas in L1544 is a = 0.19 km s -1 (Tafalla et al. 1998). 
These authors give AV = 0.22 km s _I as the typical linewidth for their observations of 
C 34 S in this region. For such a heavy molecule, turbulence is the main contributor to the 
linewidth, which allows us to estimate the mean square turbulent velocity along a typical 
direction (e.g., the line of sight) as v 2 = AV’ 2 /8In2. We easily compute that v 2 has only 
24% the value of a 2 . Assuming that it is possible to account for the ‘pressure” effects of 
such weak turbulence by adding the associated velocities in quadrature, a 2 + v 2 , we adopt 
an effective isothermal sound speed of a = 0.21 km s -1 for L1544. 

The radius of the Arecibo telescope beam at the distance of L1544 is R = 0.06 pc 
(Crutcher k Troland 2000). Ambipolar diffusion calculations by Nakano (1979), Lizano 
&: Shu (1989), Basu k Mouschovias (1994) suggest that 1 « 2 when the pivotal state is 
approached (see the summary of Li k Shu 1996). Putting together the numbers, cos z = 0.64, 
R = 0.06 pc, a = 0.21 km s l , and A = 2, we get (B\\) = 11 /zG, in excellent agreement 
with the Zeeman measurement of Crutcher k Troland (2000). These authors also deduce 
ii = 1-8 x 10 22 cm 2 from their OH measurements, whereas we compute a hydrogen column 
density within the Arecibo beam of Nh = 1.4 x 10 22 cm -2 . The slight level of disagreement is 
probably within the uncertainties in the calibration or calculation of the fractional abundance 
of OH in dark clouds (cf. Crutcher 1979, van Dishoeck k Black 1986, Flower 1990, Heiles 
et al. 1993). 

Our ability to obtain good fits of much of the observational data concerning the prestel- 
lar core L1544 with a simple analytical model should be contrasted with other, more elab- 
orate, efforts. Consider, for example, the axisymmetmc numerical simulation of Ciolek & 
Basu (2000), who were forced to assume a disk close to being edge-on (cosz % 0.3 when 
e is assumed to be 0) to reproduce the observed elongation, but who left unexplained the 
eccentric displacement of the cloud core’s center (very substantial for ellipses of eccentricity 
e % 0.54). The adoption of axisymmetric cores leads to another problem: Ciolek k Basil’s 
deprojected magnetic field is on average 3-4 times stronger than ours, values never seen di- 



rectly in Zeeman measurements of low-mass cloud cores. [See the comments of Crutcher & 
Troland (2000) concerning the need for magnetic fields in Taurus to be all nearly in the plane 
of the sky if conventional models are correct.) Natural elongation plus projection effects, as 
anticipated in the comments of Shu et al. (1999), allow us to model L1544 as a moderately 
supercritical cloud, with A ss 2, fully consistent with the theoretical expectations from am- 
bipolar diffusion calculations, and in contrast with the value A ~ 8 estimated by Crutcher 
& Troland (2000) from the measured values of B\\ and N H ■ In addition, if L1544 is a thin, 
intrinsically eccentric , disk seen moderately face-on, as implied by our model, then the ex- 
tended inward motions observed by Tafalla et al. (1998; see also Williams et al. 1999) may 
be attributable to a (relatively fast) core-amplification mechanism that gathers gas (neutral 
and ionized) dynamically but subsonically along magnetic field lines on both sides of the 
cloud toward the disk’s midplane. 

Finally, we show in Figure 2 the direction of the average magnetic field projected in 
the plane of the sky predicted by our model (thin solid line) and derived from submillimeter 
polarization observations of Ward-Thompson et al. (2000) (thin dashed line). Since we have 
assumed in our model that the major axis of iso-surface-density contours is in the plane of 
the sky, the predicted projection of the magnetic field is parallel to the cloud’s minor axis. 
The offset between the measured position angle of the magnetic field and the cloud’s minor 
axis might indicate some inclination of the cloud s major axis with respect to the plane of 
the sky. The turbulent component of the magnetic field, not included in our model, may 
also contribute to the observed deviation. 


4. Linear Perturbations of Axisymmetric Rotating SIDs 


We now consider equilibrium configurations with internal motions: 9 / 0, U(<p) * 0. 
For comparison with the analysis of Paper I, we begin with a perturbative analysis of the 
equations of the problem valid for small deviation from axisvmmetry. 

For axisymmetric disks, V m = 0 for m > 1, and therefore equations (40) and (41) reduce 
to 

S = So = 1 , V — Vo = — In 2. (53) 

Iso-surface-density contours are now circles. Substitution of these values into equation (33) 
and (34) yields B = D and U = 0. The dynamics of centrifugal balance is contained in the 
relationship (30) among the various constants of the problem: 


I< = 


©a 2 

2ireG 


(1 +D ! ), 


(54} 


the same as equation (9) of Paper I. These axisymne ■ ;c SIDs are uniquely determined, in 
an irreducible sense, by a freely specifiable value of D. bysically, we are also free, of course, 
to choose different scalings via a and A, with the latt* determining e and 0.) 

Consider now small departures from these axisyi. metric states characterized by a basic 
yV/-fold symmetry, with M = 1,2,3, . . .. Equations ( 1 and (41) give 


S((f) — 1 — MV\t cos 

••-V), 

(55) 

V{ip) = — In 2 + Vm a 

•M- 

(56) 

Equation (55) shows that for small deviations from ax. 
are limagons of Pascal (Diirer 1525). 

mrnetric iso-surface-density contours 

As required by equation (32), U((p) must be exp 

:ded as a sine series, 


U(<p) = Um sin(Ai 


(57) 

To linear order, B = D as in the axisymmetric case. 



Substitution of the relations (55)-(57) into equal ... 
yields, after subtraction of the axisymmetric relatio: 

as (34) and (32) of the governing set 
ad linearizing, 

A/ 2 (l - D 2 )V M - M{ 1 + d 2 y 

— DU m = 0, 

(58) 

II 

+ 

Q 

1 


(59) 

Solutions are possible for arbitrary (infinitesimal) va!u. 

•s of V M provided 


U M = 2D\\,. 


(60) 

and 

M( 1 + D 2 ) - M 2 { 1 - D- i 

= 2D 2 . 

(61) 

Equation (61) is equivalent to equation (25) of Paper i 
any rotation rate D (including D = 0). For M >1, we 

I and can be satisfied by M = 
require special values of D: 

1 for 

£> 2 = M for M = 

M + 2 

2.3,4,... 

(62) 

Notice the result that the required D 2 — ► 1 as M — > ex; 



For any given Z), different values of V M <C 1 generate a continuum of linearized solutions. 
Without loss of generality, we can assume V AI > 0, as the transformation V M — > — V Af is 
equivalent to a rotation of the equilibrium configuration by an angle n/M. (see discussion 



at the end of §2.2). To lowest order, the two components of the fluid velocity as given by 
equation (29) satisfy 


0'/ 2 a 


= 2DVm sin(Mip), 


and 


u 


v> _ 


= D[ 1 + V M cos {M<p)\. 


(63) 


(64) 


0 1 / 2 a 

Therefore, for infinitesimal values of V M the flow describes a locus in the velocity plane 
(u W) u, fi ) which is an ellipse of axial ratio 2 centered on (O,0^ 2 aZ)): 


+ (-^~ 
V2 0‘/ 2 J V0 1 / 2 a 


D ) 2 = D 2 V 2 


M- 


(65) 


Notice that the axial ratios are a factor of \/2 larger than the kinematic epicycles a colli- 
sionless body would generate upon being disturbed from a circular orbit in a disk that has 
a flat rotation curve (e.g., Binney & Tremaine 1987); the extra factor of \/2 (and a non- 
precessing pattern with M lobes) arises for a fluid disk because of the coherence enforced by 
the collective self-gravity of the perturbations. 


As V\i is increased, the flow must eventually try to cross the magnetosonic point, 
Up = G l / 2 aD, which is a singular point of equation (34). This transition cannot be followed 
without the introduction of shocks (see the analogous phenomena of spiral galactic shocks 
treated by Shu, Milione, & Roberts 1973). In the present context, smooth-flow solutions 
are possible only if Up < aO l ^ 2 (entirely submagnetosonic flow, for D < 1) or u.p > a© 
(entirely supermagnetosonic flow, for D > 1). When D is close to 1, either slightly smaller or 
larger, the azimuthal velocity in the SID is very close to magnetosonic already in the axisym- 
metric case. Thus, the magnetosonic point is reached when deviations from axisvmmetry 
are small, and the results of the linear analysis developed above can be applied. Equation 
(64) then gives the critical value of the coefficient V\/, in the linear regime, at which the flow 
tries to cross the magnetosonic point, 



with the plus (minus) sign valid for D > 1 (D < 1). 


5. Fully Nonlinear Models with Internal Motions 
5.1. Numerical Method 

In the general case, we solve the set of governing equations by iteration. For a given 
iterate when S(<p) is known, we may regard equation (23) as an integral for V(ip). Similarly 


- 21 


equations (32) and (28) constitute an ODE plus its starting condition for U{p). For general 
p, equation (34) may now be solved as a first order ODE with the boundary condition 
equation (16) to obtain a new iterate for S(</?). The procedure actually adopted substitutes 
a Fourier transform for a direct integration of equation (23), as described in §2.3. 

(A) Fix the value of D that one wants to study. Suppose we want to study a configuration 
with a basic A/-fold symmetry, with M = 1,2,3,.... Then we would begin with an initial 
guess for the Fourier coefficients We then compute 

OO 

V(<p) = - In 2 + V m cos(rnMp), (67) 

m— 1 

and 

OO 

S(p) = 1 - M mV m cos(mM<p). (68) 

TTl— 1 


(B) Compute the resulting value of B from equation (33). Since the cycle need be taken 
only over 2ir/M in p, we have 


MD r**/u dp 
'h r Jo S(p) 


(69) 


Integrate equation (32) for U subject to the starting condition (28). Since S has been forced 
to be a cosine series, U is then automatically a sine series, i.e., we should automatically find 
U(p) to be M-periodic, with U{2njM\ — 0. 

(C) With D fixed, and with B , V(p), S(p), and U(p), known in the form of the current 
iterates, solve equation (34) as a first order ODE for S(<p), subject to the normalization 
condition (16). With this new iterate for S((p) compute the Fourier coefficients 


\ r2nj\i 

V m = -—I S(p) cos(mMp) dp for m = l,2, ... (70) 

Compare these coefficients with those from the previous iterate. If they are insufficiently 
precise, go back to step (A), after introducing, if necessary, a relaxation parameter to smooth 
between successive iterates for V m . 


5.2. Numerical Results 

Results from our numerical integrations are illustrated in Figures 3-10. It is convenient 
to define a plane (D 2 , S,\(), where S m = —M\' m is the first coefficient in the Fourier expansion 
of the function S(p), and can be considered an indicative measure of deviations from axial 



symmetry. Figure 3 shows the regions in the ( D 2 , Si) plane occupied by M — 1 models with 
entirely submagnetosonic or entirely supermagnetosonic flow. At the upper limit of these 
two regions the flow attempts a magnetosonic transition at perisys (closest to the system 
center) in the former case and at aposys (farthest from the system center) in the latter case, 
as computed numerically with the method described in §5.1. The long-dashed line shows 
the same magnetosonic limit as given by equation (66) in the linear approximation Si 1. 
Notice that for D = 0 the results of §3 show that Sf' 1 = V™ 1 = 2. Tick marks denote the 
values of D 2 , as predicted by the linear analysis of Paper I and §4, where bifurcations occur 
with A/-fold symmetry (M > 2) from the axisymmetric sequence of SIDs that lie along the 
short dashed line. 

Figure 4 shows submagnetosonic M = 1 states for the case D 2 = 0.1 as Si progresses 
from the axisymmetric limit (Si = 0) to just before the magnetosonic transition (Si = 1.39). 
Notice that flow velocities are largest at perisys because of the tendency to conserve specific 
angular momentum (not exact because the self-consistent gravitational field is nonaxisym- 
metric). As a consequence, the magnetosonic transition, when it arrives, is made at the 
minimum of the gravitational potential, as seen by a fluid element, when the base flow is 
submagnetosonic. Notice also that the iso-surface-density contours are quasi-elliptical with 
foci at the center of the system and with the major axes lying in the same direction as the 
elongation of the streamlines formed by connecting the flow arrows. 

Figure 5 shows supermagnetosonic M = 1 states for the case D 2 = 4 as Si progresses 
from the axisymmetric limit (S i = 0) to just before the magnetosonic transition (Si = 1.08). 
Notice that flow velocities are smallest at aposys, again because of the (inexact) tendency to 
conserve specific angular momentum. As a consequence, the magnetosonic transition, when 
it arrives, is made at the maximum of the gravitational potential, as seen by a fluid element, 
when the base flow is supermagnetosonic. Notice also that the iso-surface-density contours 
are now elongated in the opposite sense to streamlines made by connecting the flow arrows. 

We can explain the last difference between the submagnetosonic and supermagnetosonic 
cases (compare Figs. 3 and 4) by analogy with a forced harmonic oscillator, whose response 
is in phase or out of phase with the external sinusoidal forcing depending on whether the 
forcing frequency is lower or higher than the natural frequency. A similar effect evidently dis- 
tinguishes the ability of fluid elements to respond in or out of phase to the nonaxisymmetrie 
forcing of the collective gravitational potential depending on whether the flow occurs at sub- 
magnetosonic or supermagnetosonic speeds relative to the pattern speed (zero in the present 
case). This distinction could be developed as a powerful diagnostic of physical conditions in 
flattened cloud cores and massive protostellar disks, if both turn out to have lopsided shapes, 
because the former can generally be expected to have submagnetosonic rotation speeds; the 


latter, supermagnetosonic speeds. 

Figure 6 shows additional examples of entirely s\ aagnetosonic flow (for D 2 = 0.5) and 
entirely supermagnetosonic flow (for D 2 = 1.5) for --- 1 SIDs, but now in the (a B ,u y ) 
plane. Models are computed with different values of, <v the numerical method described in 
§5.1. For comparison, the corresponding flow solutioi >btained with the linear analysis of §4 
are also shown. Notice that the forced epicyclic moti< >v the nonaxisymmetric gravitational 
field about the gyrocenter marked with a cross (c-. vsponding to circular motion of the 
axisymmetric model with the same value of D 2 ), ;> roaches the magnetosonic transition 
(horizontal dashed line) in both cases along a t,ang< : in the velocity-velocity plane. This 
behavior is peculiar to M = 1 SIDs, and constitute ;< topic to which we will return after 
discussing the M > 1 cases. 

Figure 7 shows the locus in the D 2 —\Sm\ plan' ‘4 sequences of equilibria with given 
A/- fold symmetry, ranging from axisymmetric mode! dashed line) to the points where the 
submagnetosonic flow acquires a magnetosonic tran.-vion (circles). We remind the reader 
that, unlike the A/ = 1 case, bifurcation of M > 1 nonces from the axisymmetric state 
occurs at discrete rather a continuum of values of /> . given bv D 2 = M/(M 4- 2). Thus, 
M = 2,3,4,... sequences always begin submagne: -onically, D 2 < 1, at Sm = 0, and 
terminate with a magnetosonic transition (circles) b< • .* the nonlinearitv parameter Sm can 
acquire very large values. 

Figure 8 shows iso-surface-density contours and -docity vectors for M = 2 equilibria 
ranging from the axisymmetric limit ( S 2 = 0) to juM before the magnetosonic transition 
(S 2 = 0.229). Notice the transformation from oval distortions at small S 2 (e.g., S 2 = 0.1) 
to dumbells at large S 2 (e.g., S 2 = 0.2). The latter shapes terminate at the magnetosonic 
transition (S 2 = 0.229), where the pinched neck ot trie dumbell develops a cusp and the 
streamlines are trying to change from circulation around a single center of attraction to 
circulation around what looks increasingly like two centers of attraction. 

Figure 9 shows iso-surface-density contours and streamlines for models with A/-fold 
symmetry. A/ = 2, 3, 4 and 5, near the endpoints of the sequences shown in Figure 7. Finally, 
Figure 10 shows the velocity-velocity plots for the same four models. The solutions with 
A/ > 1 in Figure 10 differ from those with A/ — 1 in Figure 4 in that the magnetosonic 
transition for M > 1 are made via the development of a cusp in both the iso-surface-density 
and velocity-velocity plots. We noted earlier that the magnetosnic transition is made for 
A/ = 1 configurations with the u w — u 9 locus becoming tangent to the critical curve. 



5.3. Interpretation as Onset of Shocks 


For gas flow in spiral galaxies, Shu et al. (1973) identified cusp formation in the velocity- 
velocity plane, as the onset of a shockwave with infinitesimal jumps, and we adopt a similar 
interpretation here. For trans-magnetosonic flow beyond the cusp solution (not shown in 
Figure 10 but see Shu et al. 1973), a smooth transition from submagnetosonic speeds to 
supermagnetosonic speeds is possible as the gas swings toward its closest approach to the 
center, but a smooth deceleration from supermagnetosonic speeds back to submagnetosonic 
speeds is not possible as this gas climbs outwards and catches up with slower moving material 
ahead of it. The transition to slower speeds is made instead via a sudden jump (a shockwave 
of finite strength). The shock jump introduces irreversibility to the flow pattern. Prior 
to the appearance of the shockwave, the flow can equally occur in the reverse direction as 
in the forward direction, and the streamlines close on themselves. After the appearance 
of a shockwave, time reversal is no longer possible, and the streamlines no longer close 
(see, e.g., the discussions of Kalnajs 1973 and Roberts & Shu 1973). Instead, angular 
momentum is removed from the gas (via gravitational torques when the patterns of density 
and gravitational potential show phase lags) and transferred outward in the disk, causing 
individual streamlines to spiral toward the center and increasing the central concentration 
of mass. The problem then becomes intrinsically time-dependent and cannot be followed by 
the steadv-How formulation given in the present paper. 

We are uncertain why the magnetosonic transition in the case M = 1 is not made via 
cusp-formation. It may be that in this special case, sufficient gravitational deceleration from 
supermagnetosonic to submagnetosonic speeds (rather than via pressure forces) can occur 
as to allow a smooth trans-magnetosonic flow to occur in a complete circuit. Unfortunately, 
we are unable to study this unprecedented behavior by the methods of the present paper 
because the numerical errors introduced by the truncated Fourier treatment of Poisson's 
equation compromise our ability to judge true convergence in these difficult circumstances. 
In any case, it is hard to believe, even if smooth trans-magnetosonic solutions could be 
found for lopsided SIDs, that such solutions could be stable (in a time-dependent sense) to 
the creation of shockwaves by small departures from perfect 1-fold symmetry. 


5.4. Circulation and Energy 

It is interesting to ask whether the nonaxisymmetric bifurcation sequences studied in 
this paper represent merely adjacent equilibria, or also possible evolutionary tracks that 
might be accessed by secular evolution of a single system. To help answer this question, it is 
useful to compute the variation of four quantities along any sequence. The first quantity is 


the ratio C = C/ M of the circulation C associated with a streamline to the mass M that it 
encloses. For scale-free equilibria, the value of C is independent of the spatial location of the 
streamline used to perform the calculation. The second, third, and fourth quantities are the 
ratios T = T/M, P = V/M, and W = W/M, respectively, of the kinetic energy T, pressure 
work integral P, and gravitational work integral W contained interior to any streamline, to 
the enclosed mass M. The quantity C is interesting because Kelvin's circulation theorem 
(e.g., Shu 1992) combined with the equation of continuity states that C is conserved in 
any time-dependent evolution of an ideal barotropic fluid. The quantities T, P, and W are 
interesting because they must satisfy the following scalar virial theorem (per unit mass): 

2T + QP + eW = 0. (7 1 \ 


Let m - m 0 (<p) define a streamline in the plane of the disk. The condition T = const 
in equation (24) gives immediately 


&o(<p) cx e W(v)/D , 


(72) 


where the value of the proportionality constant is irrelevant for what follows [the reader 
should not confuse the function W(<p) with If = W/M}. The mass and kinetic energy- 
contained interior to this streamline are 


M.rj 

1 f 2 * r~ o{y) 


’2ir pro q(^) 


7 dzu d<p . 
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whereas the circulation and pressure and gravitational w-ork integrals associated with this 
streamline are 




v = -f fx vnd 2 l = j 

= - j j Zx VVd'x = - f dv j/° M 


dm 

■ruo(ip) Q\? 


(75) 

(76) 

(77) 

Notice that the quantity V equals twice the thermal energy minus a surface term only if we 
perform an integration by parts, which we do not do here (cf. §3.2 in Paper I). 

If we introduce the nondimensional variables defined in §2.2, these expressions become 

T=^KQa 2 I 2 . 


M = KI { . 


(78) 



(79) 
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C ~ D 2 ' 


V = a 2 A'/,, 


W = -2nGK 2 1\, 


where we have used equation (22) to evaluate tzdV/dzu as 2nGK, and where 
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Multiplying equation (32) by tu 0 ((p) defined by equation (72) and integrating over a complete 
cycle, we obtain 

Y = DB. (81) 

A 
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and 
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where we have used equation (30) to eliminate K. With the expressions (83), the scalar 
virial theorem (71) is satisfied identically. 

Since M = 1 equilibria exist as a densely populated set of points in the D 2 -|Si| plane, 
it is clear that we can choose many sequences for them that have constant values for C. For 
fixed A (field freezing) and a (isothermal systems), C is constant along curves of constant 
B/( 1 + DB) = Do /(l + Dq), where D 0 is the axisvmmetric value of D. Thus, on such a 

sequence, 

D 0 D 


BD = 


1 + Dq{Dq — D) 


(84) 


The dotted curves in Figure 3 show such loci for two representative sequences in the D 2 - 
|Si| plane: one submagnetosonic, the other supermagnetosonic. At the beginning and end of 
the supermagnetosonic sequence displayed in Figure 3, Dq — 1.50 and D 2 — 1.84. Hence, BD 
varies from 1.50 at the beginning to 1.98 at the end, and —W oc (1 + BD) therefore increases 
by a about 29% from beginning to end. In other words, rapidly rotating, self-gravitating 
SIDs with diplaced centers are more gravitationally bound than their axisymmetric coun- 
terparts. In the presence of dissipative agents that lower the energy while preserving the 
circulation, such disks will secularly tend toward greater asymmetric elongation (see Fig. 5). 
More gravitational energy is released when distorted streamlines bring matter closer to the 
center than is expended when the same streamlines take the matter farther from the center, 
conserving circulation. This exciting result deserves further exploration both theoretically 
and observationally for systems other than the full singular isothermal disk. 


At the beginning and end of the submagnetosonic sequence displayed in Figure 3, Dq = 
0.60 and D 2 = 0.35. Hence, BD varies from 0.60 at the beginning to 0.40 at the end, and 
-W oc (1 + BD) therefore decreases by about 12% from beginning to end. This variation 
is not very much considering how fast this sequence rotates relative to realistic cloud cores. 
Nevertheless, the formal decrease of — W as one leaves the axisymmetric state implies that 
submagnetosonic systems require some input of energy to make them less round. Exceptions 
are sequences that branch from smaller values of Z?q, which have smaller variations of — W . 
In particular, long spindles have no binding energy disadvantage whatsoever relative to 
axisymmetric disks for the nonrotating sequence shown in Figure 1, because here —IF oc 
(1 + BD) = 1, a constant. In this regard, it may be significant that observed cores that are 
significantly lopsided (see Fig. 2) typically rotate quite slowly. 

The story is more ambiguous for M > 1 equilibria. Here, for given A/, the stationary 
states occupy one-dimensional curves in the D 2 -\Sm\ plane; therefore, we have no control 
over how C and —W vary along any sequence. Plotted in Figure 11 are the values of C 
and —IT as we vary Sm along the sequences for M = 2, 3, 4, 5. Amazingly, the normalized 
circulation C is nearly, but not exactly, constant on each sequence, varying by no more than 
1% in all cases. Given the small values of Sm for which solutions exist and the relatively 
small variation of D along each sequence, this result is not surprising, because B and DB 
differ from their values for axisymmetric SIDs by terms 0{S \ t ) . Although in principle secular 
evolution along any M > 1 sequence would require a slight redistribution of circulation with 
mass, the amount required is truly slight, and one could imagine that mechanisms might 
exist that can effect a slow transformation along the sequence toward more nonaxisymmetric 
states. In principle, such evolution would seem to favor the formation of M = 2, 3, 4, 5. 

■ * * buds, depending on the rate of rotation present in the underiving flow. However, before 2, 
3, 4, 5, . . . independently orbiting bodies can form by such a "fission” process, this sequence 
of events would terminate in shockwaves, and the resultant transfer of angular momentum 
(or circulation) outward and mass inward would stabilize the system against actual successful 
fission. 

In practice, for gaseous systems, a more practical difficulty mitigates against even be- 
ginning the secular paths of evolution described in the previous paragraphs for the submag- 
netosonic cases. The nonaxisymmetric SIDs with M = 2, 3, 4, 5 depicted in Figures 8 and 9 
are all rotating too slowly to be stable against “inside-out” collapse of the type studied for 
their axisymmetric counterpart by Li & Shu (1997). This dynamical instability would for- 
mallv overwhelm any secular evolution along the lines described above. (Supermagnetosonic 
M _ 1 configurations rotate quickly enough to be stable against “inside-out” dynamical col- 
lapse, and a secular transformation to the more elongated and eccentrically displaced states 
of Fig. 5 are realistic theoretical possibilities.) We plan to study the dynamical collapse 


and fragmentation properties of nonaxisvmmetric, submagnetosonie SIDs with general M- 
fold symmetry in a future paper. In another treatment, we shall also discuss the question 
whether configurations with strict M > 1 symmetry are formally (secularly) unstable also 
to perturbations of A / = 1 periodicity (i.e., to additional lopsided bifurcations). But, for 
the present, we merely remark that the practical attainment of any of the nonaxisymmet- 
ric pivotal states depicted, say, in Figure 8 probably occurs, not along a sequence where 
each member has already achieved a (nearly) singular value of surface density at the origin 
w = 0, but along a line of evolution (perhaps by ambipolar diffusion) where the growing 
central concentration of matter occurs without the a priori assumption of axial symmetry 
(e.g.. nonaxisvmmetric generalizations of the calculations of Basu & Mouschovias 1994). 


6. Summary and Discusssion 

In this paper we have shown that prestellar molecular cloud cores modeled in their piv- 
otal state just before the onset of gravitational collapse (protostar formation and envelope 
infall) as magnetized singular isothermal disks need not be axisymmetric. The most impres- 
sive distortions are those that make slowly rotating circular cloud cores lopsided (M = 1 
asvmmetrv). Although slowly rotating, lopsided cloud cores have a slight disadvantage rel- 
ative to their axisymmetric counterparts from an energetic point of view, such elliptical 
configurations do seem to appear in nature (see Fig. 2). 

More intriguingly, elongated, eccentrically displaced, supermagnetosonicallv rotating 
SIDs (that are stable to overall graviational collapse) are preferred to their axisymmetric 
counterparts if the excess binding energy of the latter can be radiated away without chang- 
ing the circulation of the streamlines. If the mass of the circumstellar disk of a very young 
protostar is a large fraction of the mass of the system, it might be possible to find such 
M = 1 distortions of actual objects by future MMA observations. If such disks have (per- 
turbed) flat rotation curves, we predict (see Fig. 5) that the mm-wave isotphotes should be 
elongated perpendicular to an eccentrically displaced central star and also 'perpendicular to 
the eccentric shape of the streamlines (as might be deducible from isovelocity plots common 
for investigations of spiral and barred galaxies). 

Bifurcations into sequences with M = 2, 3, 4, 5, and higher symmetry require rotation 
rates considerably larger (>0.7 times the magnetosonic speeds) than is typically measured 
for observed molecular cloud cores (e.g., Goodman et al. 1993). Although seemingly more 
promising for binary and multiple star-formation, the models with M = 2, 3, 4, 5, ... 
symmetries all terminate in shockwaves before their separate lobes can succeed in forming 
anything that resembles separate bodies (see Fig. 8). For these configurations to exist at all, 


the basic rotation rate has to be fairly close to magnetosonic. It is then not possible for the 
nonaxial symmetry to become sufficiently pronounced as to turn streamlines that circulate 
around a single center to streamlines that circulate around multiple centers (as is needed 
to form multiple stars), without the distortions causing supermagnetosonically flowing gas 
to slam into submagnetosonically flowing gas. The resultant shockwaves then increase the 
central concentration in such a fashion as to suppress the tendency toward fission. 

VVe have managed to gain the above understanding semi-analytically only because of the 
mathematical simplicity of isopedically magnetized SIDs. The same understanding probably 
underlies similar findings from numerical simulations of the fission process that inevitably end 
with the creation of shockwaves before the actual production of two or more separately grav- 
itating bodies (Tohline 2000, personal communication). This negative result, combined with 
the analysis of the spiral instabilities that afflict the more rapidly rotating, self-gravitating, 
disks into which more slowly rotating, cloud cores collapse (also modeled here as SIDs), is 
cause for pessimism that a successful mechanism of binary and multiple star-formation can 
be found by either the fission or the fragmentation process acting in the aftermath of the 
gravitational collapse of marginally supercritical clouds dunng the stages when field freezing 
provides a good dynamical assumption. 

It might be argued that our analysis also assumed smooth starting conditions, and that 
theiefoie, turbulence might be the more important missing ingredient. However, the low- 
mass cloud cores in the Taurus molecular cloud that gives rise to many binaries and multiple- 
star systems composed of sunlike stars are notoriously quiet, with turbulent velocities that 
are only a fraction of the thermal sound speed (e.g., Fuller Sz Myers 1992). Such levels of 
turbulence are well below those that appear necessary to induce “turbulent fragmentation” 
in the numerical simulations of Klein et al. (2000). Interstellar turbulence is undoubtedly 
an important process at the larger scales that characterize the fractal structures of giant 
molecular clouds (see, e.g., Allen & Shu 2000), but it probably plays only a relatively minor 
role in the simplest case of isolated or distributed star-formation that we see in clouds like 
those in the Taurus-Auriga region, which has, as we mentioned earlier, more than its share 
of cosmic binaries. 


In contrast, we know that the dimensionless mass-to-flux ratio A had to increase from 
values typically ~ 2 in cloud cores to values in excess of 5000 in formed stars (Li & Shu 1997). 
Massive loss of magnetic flux must have occurred at some stage of the gravitational collapse 
of molecular cloud cores to form stars. Moreover, this loss must take place at some point 
at a dynamical rate, or even faster, since the collapse process from pivotal molecular cloud 
cores is itself dynamical. It is believed that dynamical loss of magnetic fields from cosmic 
gases occurs only when the volume density exceeds ~ 10“ H 2 molecules cm' 3 (e.g., Nakano 


k Lmebayashi 198Ga,b; Desch k Mouschovias 2000). It might be thought that cloud cores 
have to collapse to fairly small linear dimensions before the volume density reaches such high 
values, and therefore, that only close binaries can be explained by such a process, but not 
wide binaries (McKee 2000, personal communication). However, this impression is gained 
by experience with axisymmetric collapse. Once the restrictive assumption of perfect axial 
symmetry is removed, we gain the possibility that some dimensions may shrink faster than 
others (e.g., Lin, Mestel, k Shu 1965), and densities as high as 10“ cm -3 might be reached 
while only one or two dimensions are relatively small, and while the third is still large enough 
to accomodate the (generally eccentric) orbits of wide binaries. 

We close with the following analogies. The basic problem with trapped magnetic fields 
is that they compress like relativistic gases (i.e., their stresses accumulate as the 4/3 pow- 
er increase of the density in 3-D compression). Such gases have critical masses [e.g., the 
Chandrasekhar limit in the theory of white dwarfs, or the magnetic critical mass of equa- 
tion (2)] which prevent their self-gravitating collections from suffering indefinite compression, 
no matter how high is the surface pressure, if the object masses lie below the critical values. 
Moreover, while marginally supercritical objects might collapse to more compact objects 
(e.g., white dwarfs into neutron stars, or cloud cores into stars), a single such object cannot 
be expected to naturally fragment into multiple bodies (e.g., a single white dwarf with mass 
slightly bigger than the Chandrasekhar limit into a pair of neutron stars). 

In order for fragmentation to occur, it might be necessary for the fluid to decouple 
rapidly from its source of relativistic stress. For example, the universe as a whole always 
has many thermal Jeans masses. Yet in conventional big-bang theory, this attribute did 
not do the universe any good in the problem of making gravitationally bound subunits, 
as long as the universe was tightly coupled to a relativistic (photon) field. Only after the 
matter field had decoupled from the radiation field in the recombination era, did the many 
fluctuations above the Jeans scale have a chance to produce gravitational “fragments." It is 
our contention that this second analogy points toward where one should search for a viable 
theory of the origin of binary and multiple stars from the gravitational collapse of magnetized 
molecular cloud cores. 
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by a grant from the National Science Foundation and in part by the NASA Astrophvsical 
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A. Appendix 

Given the definition 

^ m ~ 9^ — cos cos { T12i p) dp = — J ln(l — cos p) cos(rrup) dp, 

where m is an integer (positive or negative), we show that 

^ f — In 2 if m = 0, 

m ~ \ — l/|m| if |mi > 1. 

(1) For m = 0, 

L ° = 2 ^ / ln ^ “ C0S ^ d{p ~ ~ { ,n ( 1 ~ C0St P) di P 
Successive transformations p -> -p and p -4 p + - demonstrate 

If 0 .. 1 r' 
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1 / 1 /*T 

= -/ ln(l-cos p)dp~-j ln(l + cos p) dp. 

” J 7 r i(J 


(A3) 


(A4) 


which is just a statement that cosy? is odd relative to the midpoint of the interval (0, 7 r). If 
we add equations (A3) and (A4), we get 

2 ^° = 7 t Jo ln ( 1 ~ cos2 ’-p) dp = - j ln(sin 2 p) dp. (A5) 

Rewrite sin 2 ^ = (1/2)(1 - cos \) where x = 2 p: then 

1 f 2 - 

2L 0 = — In 2 4- — \n(l - cos x) d\. (A 6) 

But the last integral is another expression for L 0 (see eq. [A3]); thus, L 0 = - In 2. (QED) 

(2) Integrate equation (Al) once by parts to obtain 

T _ 1 [+* sin [mp) ,;in pdp 

m ~ 2tt L ~T^7o~' ( A7 ) 

Multiply and divide the integrand by 1 + cos p and write 1 - cos 2 p as sin 2 p. Thereby obtain 

m ^m = — (I, n + J m ), (A8) 
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Fig. 1. Iso-surface-density contours for static SIDs are perfect ellipses with eccentricity e. 
For e = 0 the SID is an axisymmetric disk with surface density oc zu~ l , for e = 1 the SID 
degenerates in a semi-infinite filament of constant linear mass density. 
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Fig. 2.— Iso-surface-brightness contours ( thick solid lines) from a theoretically computed, 
lopsided, magnetized, self-gravitating figure of equilibrium compared with isophotal mea- 
surements of Ward-Thompson et al. (1999) of the submillimeter emission from heated dust 
grains in L1544. The short solid line and dashed line show the directions of predicted and 
measured field inferred from submillimeter-wave polarization observations (Ward-Thompson 

et al. 2000). 
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Fig. 3. — The loci (solid curves) in the ( D 2 , |Si|) plane of critical flow solutions for M = 1 
approaching magnetosonic speed starting from entirely submagnetosonic ( D < 1) or entire- 
ly supermagnetosonic flow ( D > 1). The horizontal short-dashed line shows the locus of 
axisymmetric models. The long-dashed curves shows the limit of magnetosonic models as 
given by equation (66) in the linear approximation 5] 1. The dotted curves show a sub- 

magnetosonic and a supermagnetosonic sequence, each of which maintains a constant ratio 
C of circulation C to enclosed mass M (see §5.4). iVkinarks pointing downward from the 
horizontal dashed line denote the values of D 2 where distortions with M-fold symmetry can 
occur, with M > 1, as predicted by the linear analysis of Paper I and §4. 
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Fig. 4 — Sequence of equilibria with M = 1 and D 2 = 0.1 for S x = 0 (axisymmetric case). 
0.3, 0.6 and 1.39 (magnetosonic flow). Thick lines correspond to iso-surface-density contours 
logarithmically spaced. Streamlines are outlined by vectors of length proportional to the 
modulus of velocity, drawn to scale in the four panels. Equation (29) shows that the flow 
velocities depend, in magnitude and direction, only on the azimuthal angle p and not on 
the distance zd (for given p) from the center of the mass distribution. In particular, the gas 
velocities at perisys (or aposys) are all the same independent of the size of the streamline. 
In the last panel, the thin dotted line indicates the locus of magnetosonic velocity reached at 
perisys. 
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Fig. 6. — Two examples of entirely submagnetosonic flow ( D 2 = 0.5) and entirely super- 
magnetosonic flow ( D 2 = 1.5) in the velocity- velocity plane, for M = 1. In each panel, the 
axisymmetric solutions is marked by a cross. The two inner solid curves are obtained with 
Si = 0.1 and 0.2, whereas the dotted curves show the corresponding linearized solutions 
The outermost solid curve., obtained with Si = 0.3235 for D ‘ = O.o and Si = 0.2326 foi 
D 2 = 1.5, shows the approach to magnetosonic flow ( dashed line) defined by the condition 
|u| 2 = 0a 2 : 
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Fig. 7.— Locus in the D 2 -\S M \ plane of sequences of equilibria with given M-fold symmetry. 
The dashed line indicates the locus of axisymmetric equilibria. Tickmarks denote the values 
of D 2 where distortions with M-fold symmetry can occur, as predicted by the linear analysis 
of Paper I and §4. Circles indicate the points where the sequences terminate because of the 
occurrence of shocks. 







Fig. 9. Iso-surface-density contours ( thick lines) and streamlines ( thm lines) for models 
with Af-fold symmetry M = 2,3,4 and 5 at the point of shock formation. Iso-surface-densitv 
contours develop cusps where the azymuthal flow velocity reaches the magnetosonic value.' 
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